c --------------------- c subroutine initvrml c --------------------- c **** write header **** subroutine initvrml write(10,10) 10 format("#VRML V2.0 utf8") return end c --------------------- c subroutine line3d c --------------------- c **** draw lines **** subroutine line3d(lasl,fc,po) c c lasl: number of points on line c r: red color [0.0,1.0] c g: green color [0.0,1.0] c b: blue color [0.0,1.0] c po(3*lasl): x,y,z-values in position c dimension po(1),fc(3) c call linebe(fc) c do 40 ii=1,lasl lin=ii li1=3*(lin-1)+1 li2=3*(lin-1)+2 li3=3*lin x=po(li1) y=po(li2) z=po(li3) call linep(x,y,z,ii,lasl) 40 continue c call lineen(lasl) return end c --------------------- c subroutine linep c --------------------- c **** determine positions of lines ***** subroutine linep(x,y,z,ii,lasl) c if(ii.lt.lasl) then write(10,10) x,z,y c 10 format(f5.1,1x,f5.1,1x,f5.1,1h,) c 10 format(f5.2,1x,f5.2,1x,f5.2,1h,) 10 format(f8.2,1x,f8.2,1x,f8.2,1h,) else write(10,20) x,z,y c 20 format(f5.1,1x,f5.1,1x,f5.1) c 20 format(f5.2,1x,f5.2,1x,f5.2) 20 format(f8.2,1x,f8.2,1x,f8.2) end if return end c --------------------- c subroutine linebe c --------------------- c **** determine color of lines **** subroutine linebe(fc) c dimension fc(3) c write(10,21) 21 format("Shape {") write(10,22) 22 format("geometry IndexedLineSet {") write(10,23) 23 format("color Color {") write(10,20) fc(1),fc(2),fc(3) 20 format("color [ ",f4.1,1x,f4.1,1x,f4.1," ]") write(10,25) 25 format("}") write(10,26) 26 format("coord Coordinate {") write(10,27) 27 format("point [") return end c --------------------- c subroutine lineen c --------------------- c **** connect lines **** subroutine lineen(lasl) c integer ic(9000) c do 10 i=1,9000 ic(i)=i-1 10 continue write(10,21) 21 format("]") write(10,22) 22 format("}") write(10,23) 23 format("colorIndex [ 0 ]") write(10,24) 24 format("coordIndex [") c write(10,20) (ic(i),i=1,lasl) 20 format(10i2,90i3,900i4/(2000i5)) c write(10,25) 25 format("]") write(10,26) 26 format("colorPerVertex FALSE") write(10,27) 27 format("}") write(10,28) 28 format("}") return end c c --------------------- c subroutine dr_scale c --------------------- c **** draw outside frame **** subroutine dr_scale(fc,NGX,NGY,NGZ,D) c dimension sp(6),fc(3) c do 10 j=0,NGY,NGY do 20 k=0,NGZ,NGZ sp(1)=0 sp(2)=j*D sp(3)=k*D sp(4)=NGX*D sp(5)=j*D sp(6)=k*D call line3d(2,fc,sp) 20 continue 10 continue c do 30 j=0,NGX,NGX do 40 k=0,NGZ,NGZ sp(1)=j*D sp(2)=0 sp(3)=k*D sp(4)=j*D sp(5)=NGY*D sp(6)=k*D call line3d(2,fc,sp) 40 continue 30 continue c do 50 j=0,NGX,NGX do 60 k=0,NGY,NGY sp(1)=j*D sp(2)=k*D sp(3)=0 sp(4)=j*D sp(5)=k*D sp(6)=NGZ*D call line3d(2,fc,sp) 60 continue 50 continue return end c c --------------------- c subroutine viewpoint c --------------------- c **** change viewpoint **** subroutine viewpoint(NGX,NGY,NGZ,D) c character s s='"' c write(10,10) 0.5*NGX*D,-1.5*NGY*D,0.5*NGZ*D 10 format("\nViewpoint{position ",f8.2,1x,f8.2,1x,f8.2) write(10,20) 20 format("orientation 1 0 0 1.57 ") write(10,30) s,s 30 format("description ",a1,"entry",a1,"}") return end c c --------------------- c subroutine dr_face c --------------------- c **** draw faces **** subroutine dr_face(p,n,fc,tr) c real p(n,3) real fc(3) c write(10,10) 10 format("\nShape{") write(10,20) 20 format("appearance Appearance{") write(10,30) 30 format("material Material{") write(10,200) 200 format("ambientIntensity 0.0") write(10,40) 40 format("diffuseColor ") write(10,50) fc(1),fc(2),fc(3) 50 format(f4.1,1x,f4.1,1x,f4.1) write(10,60) tr 60 format("transparency ",f5.2) write(10,70) 70 format("} }") write(10,80) 80 format("geometry IndexedFaceSet{") write(10,90) 90 format("solid FALSE") write(10,100) 100 format(" coord Coordinate{") write(10,110) 110 format("point[") c do 1000 i=1,n call dr_facep(p,i,n) 1000 continue c write(10,120) 120 format("] }") write(10,130) 130 format("coordIndex[") do 1500 i=0,n-2 write(10,140) i 140 format(i2,",") 1500 continue write(10,160) n-1 160 format(i2) write(10,150) 150 format("] } }") return end c c --------------------- c subroutine dr_facep c --------------------- c **** determine positions of faces **** subroutine dr_facep(p,i,n) c real p(n,3) c if(i.lt.n) then write(10,10) p(i,1),p(i,2),p(i,3) 10 format(f10.2,1x,f10.2,1x,f10.2,1h,) else write(10,20) p(i,1),p(i,2),p(i,3) 20 format(f10.2,1x,f10.2,1x,f10.2) end if return end c --------------------- c subroutine mk_data c --------------------- c **** sample data-1 **** c **** distance from the origin **** subroutine mk_data(v,NX,NY,NZ,D) c real v(NX,NY,NZ) c do 10 i=1,NX do 20 j=1,NY do 30 k=1,NZ x=(i-1)*D y=(j-1)*D z=(k-1)*D v(i,j,k)=sqrt((x)*(x)+(y)*(y)+(z)*(z)) 30 continue 20 continue 10 continue return end c c --------------------- c subroutine mk_data2 c --------------------- c **** sample data-2 **** c **** distance from the center of lattice **** subroutine mk_data2(v,NX,NY,NZ,D) c real v(NX,NY,NZ) c do 10 i=1,NX do 20 j=1,NY do 30 k=1,NZ x=(i-1)*D y=(j-1)*D z=(k-1)*D v(i,j,k)=sqrt(((x-(NX-1)*D/2)**2)+((y-(NY-1)*D/2)**2) 1 +((z-(NY-1)*D/2)**2)) 30 continue 20 continue 10 continue return end c c --------------------- c subroutine mk_data3 c --------------------- c **** sample data-3 **** c **** sin function **** subroutine mk_data3(v,NX,NY,NZ,D) c real v(NX,NY,NZ) parameter (pi=3.14159) c do 10 i=1,NX do 20 j=1,NY do 30 k=1,NZ x=(i-1)*D y=(j-1)*D z=(k-1)*D v(i,j,k)=sin((x/(NX*D))*pi)+2*sin((y/(NY*D))*pi) 1 +3*sin((z/(NZ*D))*pi) 30 continue 20 continue 10 continue return end c c --------------------- c subroutine mk_data4 c --------------------- c **** sample data-4 **** c **** sin function-2 **** subroutine mk_data4(v,NX,NY,NZ,D) c real v(NX,NY,NZ) parameter (pi=3.14159) c do 10 i=1,NX do 20 j=1,NY do 30 k=1,NZ x=(i-1)*D y=(j-1)*D z=(k-1)*D v(i,j,k)=(sin((x/(NX*D))*2*pi))**2+2*sin((y/(NY*D))*pi) 1 +3*sin((z/(NZ*D))*pi) 30 continue 20 continue 10 continue return end c c --------------------- c subroutine mk_cube c --------------------- c **** separate by cube **** subroutine mk_cube(x,y,z,p,v,vt,NX,NY,NZ,D) c real p(8,3) real vt(8) real v(NX,NY,NZ) c do 100 i=0,1 p(1+(i*4),1)=x p(1+(i*4),2)=y p(1+(i*4),3)=z+i*D c p(2+(i*4),1)=x+D p(2+(i*4),2)=y p(2+(i*4),3)=z+i*D c p(3+(i*4),1)=x+D p(3+(i*4),2)=y+D p(3+(i*4),3)=z+i*D c p(4+(i*4),1)=x p(4+(i*4),2)=y+D p(4+(i*4),3)=z+i*D 100 continue c do 200 i=1,8 ix=int(p(i,1))/(int(D)) iy=int(p(i,2))/(int(D)) iz=int(p(i,3))/(int(D)) vt(i)=v(ix+1,iy+1,iz+1) 200 continue c return end c c --------------------- c subroutine rm_cube c --------------------- c **** remove the cube where face is not drawed **** subroutine rm_cube(th,vt,irm) c integer irm real vt(8) iuc=0 idc=0 c do 10 i=1,8 if(vt(i).ge.th) then iuc = iuc + 1 else if(vt(i).lt.th) then idc = idc + 1 end if 10 continue c if( (iuc.eq.8).or.(idc.eq.8) ) then irm = 1 else irm = 0 end if c return end c c --------------------- c subroutine cal_pc_sub c --------------------- c **** find existence of iso-points and positions **** subroutine cal_pc_sub(th,pt,vc,pcs,pcts) c real pt(2,3) real vc(2),pcs(3) integer pcts c if( ((vc(1).ge.th).and.(vc(2).ge.th)).or. 1 ((vc(1).lt.th).and.(vc(2).lt.th)) ) then pcts=0 do 10 i=1,3 pcs(i)=0 10 continue else pcts=1 if(vc(1).eq.th) then do 30 i=1,3 pcs(i)=pt(1,i) 30 continue return else if(vc(2).eq.th) then do 40 i=1,3 pcs(i)=pt(2,i) 40 continue return else s=vc(2)-vc(1) a=th-vc(1) b=s-a do 50 i=1,3 pcs(i)=(b/s)*pt(1,i)+(a/s)*pt(2,i) 50 continue end if end if return end c --------------------- c subroutine cal_pc c --------------------- c **** determine the number of iso-points and positions **** subroutine cal_pc(th,p,vt,pc,pct) c real pt(2,3),p(8,3),pc(12,3) real vc(2),vt(8),pcs(3) integer pct(12) integer pcts c do 30 k=0,1 do 20 j=1,3 do 10 i=1,3 pt(1,i)=p((j+(k*4)),i) pt(2,i)=p(((j+1)+(k*4)),i) 10 continue vc(1)=vt(j+(k*4)) vc(2)=vt((j+1)+(k*4)) call cal_pc_sub(th,pt,vc,pcs,pcts) pct(j+(k*4))=pcts do 50 i=1,3 pc(j+(k*4),i)=pcs(i) 50 continue 20 continue 30 continue c do 100 j=0,1 do 90 i=1,3 pt(1,i)=p((j*4)+1,i) pt(2,i)=p((j*4)+4,i) 90 continue vc(1)=vt((j*4)+1) vc(2)=vt((j*4)+4) call cal_pc_sub(th,pt,vc,pcs,pcts) pct((j*4)+4)=pcts do 80 i=1,3 pc((j*4)+4,i)=pcs(i) 80 continue 100 continue c do 200 j=1,4 do 180 i=1,3 pt(1,i)=p(j,i) pt(2,i)=p(j+4,i) 180 continue vc(1)=vt(j) vc(2)=vt(j+4) call cal_pc_sub(th,pt,vc,pcs,pcts) pct(j+8)=pcts do 190 i=1,3 pc(j+8,i)=pcs(i) 190 continue 200 continue end c c --------------------- c subroutine mk_tri c --------------------- c **** draw triangle **** subroutine mk_tri(pc,pct,fc,tr) c integer pct(12),indt(8,3) real pc(12,3),ptri(3,3),fc(3) c call mk_index_tri(indt) c do 40 i=1,8 if( (pct(indt(i,1)).eq.1).and. 1 (pct(indt(i,2)).eq.1).and.(pct(indt(i,3)).eq.1)) then do 30 k=1,3 do 20 j=1,3 ptri(k,j)=pc(indt(i,k),j) 20 continue pct(indt(i,k))=0 30 continue call dr_face(ptri,3,fc,tr) end if 40 continue return end c c --------------------- c subroutine mk_index_tri c --------------------- subroutine mk_index_tri(indt) c integer indt(8,3) c indt(1,1) = 1 indt(1,2) = 4 indt(1,3) = 9 c indt(2,1) = 2 indt(2,2) = 1 indt(2,3) = 10 c indt(3,1) = 3 indt(3,2) = 2 indt(3,3) = 11 c indt(4,1) = 4 indt(4,2) = 3 indt(4,3) = 12 c indt(5,1) = 5 indt(5,2) = 8 indt(5,3) = 9 c indt(6,1) = 6 indt(6,2) = 5 indt(6,3) = 10 c indt(7,1) = 7 indt(7,2) = 6 indt(7,3) = 11 c indt(8,1) = 8 indt(8,2) = 7 indt(8,3) = 12 c return end c c --------------------- c subroutine mk_squ c --------------------- c **** draw boxes or quadrilateral **** subroutine mk_squ(pc,pct,fc,tr) c integer pct(12),inds(3,5,4) real pc(12,3),psqu(4,3) real fc(3) c call mk_index_squ(inds) c do 40 i=1,3 do 30 j=2,5 if((pct(inds(i,j,1)).eq.1).and.(pct(inds(i,j,2)).eq.1).and. 1 (pct(inds(i,j,3)).eq.1).and.(pct(inds(i,j,4)).eq.1)) then do 20 k=1,4 do 10 m=1,3 psqu(k,m)=pc(inds(i,j,k),m) 10 continue pct(inds(i,j,k))=0 20 continue call dr_face(psqu,4,fc,tr) end if 30 continue 40 continue c do 100 i=1,3 n=0 do 90 m=1,4 if(pct(inds(i,1,m)).eq.1) n=n+1 90 continue if(n.eq.4) then do 80 k=1,4 do 70 j=1,3 psqu(k,j)=pc(inds(i,1,k),j) 70 continue pct(inds(i,1,k))=0 80 continue call dr_face(psqu,4,fc,tr) return end if 100 continue end c c --------------------- c subroutine mk_index_squ c --------------------- subroutine mk_index_squ(inds) c integer inds(3,5,4) c inds(1,1,1) = 1 inds(1,1,2) = 3 inds(1,1,3) = 7 inds(1,1,4) = 5 c inds(2,1,1) = 9 inds(2,1,2) = 10 inds(2,1,3) = 11 inds(2,1,4) = 12 c inds(3,1,1) = 2 inds(3,1,2) = 4 inds(3,1,3) = 8 inds(3,1,4) = 6 c inds(1,2,1) = 1 inds(1,2,2) = 2 inds(1,2,3) = 6 inds(1,2,4) = 5 c inds(1,3,1) = 2 inds(1,3,2) = 3 inds(1,3,3) = 7 inds(1,3,4) = 6 c inds(1,4,1) = 3 inds(1,4,2) = 4 inds(1,4,3) = 8 inds(1,4,4) = 7 c inds(1,5,1) = 4 inds(1,5,2) = 1 inds(1,5,3) = 5 inds(1,5,4) = 8 c inds(2,2,1) = 5 inds(2,2,2) = 7 inds(2,2,3) = 11 inds(2,2,4) = 10 c inds(2,3,1) = 1 inds(2,3,2) = 3 inds(2,3,3) = 11 inds(2,3,4) = 10 c inds(2,4,1) = 1 inds(2,4,2) = 3 inds(2,4,3) = 12 inds(2,4,4) = 9 c inds(2,5,1) = 5 inds(2,5,2) = 7 inds(2,5,3) = 12 inds(2,5,4) = 9 c inds(3,2,1) = 2 inds(3,2,2) = 4 inds(3,2,3) = 9 inds(3,2,4) = 10 c inds(3,3,1) = 2 inds(3,3,2) = 4 inds(3,3,3) = 12 inds(3,3,4) = 11 c inds(3,4,1) = 6 inds(3,4,2) = 8 inds(3,4,3) = 12 inds(3,4,4) = 11 c inds(3,5,1) = 6 inds(3,5,2) = 8 inds(3,5,3) = 9 inds(3,5,4) = 10 c return end c c --------------------- c subroutine count_pc c --------------------- subroutine count_pc(pct,num) c integer pct(12) c num=0 do 10 i=1,12 if(pct(i).eq.1) then num = num + 1 end if 10 continue return end c c --------------------- c subroutine mk_index_face c --------------------- subroutine mk_index_face(indf) c integer indf(6,4) c indf(1,1) = 1 indf(1,2) = 2 indf(1,3) = 3 indf(1,4) = 4 c indf(2,1) = 5 indf(2,2) = 6 indf(2,3) = 7 indf(2,4) = 8 c indf(3,1) = 1 indf(3,2) = 5 indf(3,3) = 9 indf(3,4) = 10 c indf(4,1) = 2 indf(4,2) = 6 indf(4,3) = 10 indf(4,4) = 11 c indf(5,1) = 3 indf(5,2) = 7 indf(5,3) = 11 indf(5,4) = 12 c indf(6,1) = 4 indf(6,2) = 8 indf(6,3) = 9 indf(6,4) = 12 c return end c c --------------------- c subroutine search_f c --------------------- subroutine search_f(indf,self,next,face,ptr) c integer indf(6,4) integer next,face,ptr,self c do 20 i=1,6 if(i.ne.self) then do 10 j=1,4 if(next.eq.indf(i,j)) then face=i ptr=j end if 10 continue end if 20 continue c return end c c --------------------- c subroutine search_p c --------------------- subroutine search_p(pc,pct,indf,face,ptr,pft,next) c real pc(12,3),pft(3) integer pct(12),indf(6,4) integer face,ptr,next c do 20 i=1,4 if((i.ne.ptr).and.(pct(indf(face,i)).eq.1)) then do 10 j=1,3 pft(j)=pc(indf(face,i),j) 10 continue next=indf(face,i) end if 20 continue c return end c --------------------- c subroutine mk_pent c --------------------- c **** draw pentagon **** subroutine mk_pent(pc,pct,fc,tr) c integer face,ptr,self,next real pc(12,3),pf(5,3) real fc(3),pft(3) integer pct(12),indf(6,3) c call count_pc(pct,n) if(n.ne.5) then return end if c call mk_index_face(indf) c k=1 self=1 do 10 i=1,4 if(pct(indf(self,i)).eq.1) then do 20 j=1,3 pf(k,j)=pc(indf(self,i),j) 20 continue next=indf(self,i) k=k+1 end if 10 continue c if(k.eq.1) then self=2 do 30 i=1,4 if(pct(indf(self,i)).eq.1) then do 40 j=1,3 pf(k,j) = pc(indf(self,i),j) 40 continue next=indf(self,i) k=k+1 end if 30 continue end if c do 50 k=3,n call search_f(indf,self,next,face,ptr) call search_p(pc,pct,indf,face,ptr,pft,next) self=face do 60 j=1,3 pf(k,j)=pft(j) 60 continue 50 continue c call dr_face(pf,5,fc,tr) c return end c c --------------------- c subroutine mk_hexa c --------------------- c **** draw hexagon **** subroutine mk_hexa(pc,pct,fc,tr) c integer face,ptr,self,next real pc(12,3),ph(6,3) real fc(3),pft(3) integer pct(12),indf(6,4) c call count_pc(pct,n) if(n.ne.6) then return end if c call mk_index_face(indf) c k=1 self=1 do 10 i=1,4 if(pct(indf(self,i)).eq.1) then do 20 j=1,3 ph(k,j)=pc(indf(self,i),j) 20 continue next=indf(self,i) k=k+1 end if 10 continue c if(k.eq.1) then self=2 do 30 i=1,4 if(pct(indf(self,i)).eq.1) then do 40 j=1,3 ph(k,j) = pc(indf(self,i),j) 40 continue next=indf(self,i) k=k+1 end if 30 continue end if c do 50 k=3,n call search_f(indf,self,next,face,ptr) call search_p(pc,pct,indf,face,ptr,pft,next) self=face do 60 j=1,3 ph(k,j)=pft(j) 60 continue 50 continue c call dr_face(ph,n,fc,tr) c return end c c --------------------- c subroutine mk_equ_face_sub c --------------------- subroutine mk_equ_face_sub(th,p,vt,fc,tr) c integer pct(12) real pc(12,3),p(8,3) real fc(3),vt(8) c call cal_pc(th,p,vt,pc,pct) c call mk_tri(pc,pct,fc,tr) c call mk_squ(pc,pct,fc,tr) call mk_pent(pc,pct,fc,tr) call mk_hexa(pc,pct,fc,tr) c return end c c --------------------- c subroutine mk_equ_face c --------------------- c **** draw isosurface **** subroutine mk_equ_face(th,v,fc,tr,NX,NY,NZ,D) c integer irm real v(NX,NX,NX) real vt(8),fc(3),p(8,3) c do 30 x=0,(NX-2)*D,D do 20 y=0,(NY-2)*D,D do 10 z=0,(NZ-2)*D,D call mk_cube(x,y,z,p,v,vt,NX,NY,NZ,D) call rm_cube(th,vt,irm) if(irm.eq.0) then call mk_equ_face_sub(th,p,vt,fc,tr) end if 10 continue 20 continue 30 continue c return end